An exact pressure evolution equation for the incompressible Navier-Stokes equations 
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In this paper the issue of the determination of the fluid pressure in incompressible fluids is ad- 
dressed, with particular reference to the search of algorithms which permit to advance in time the 
fluid pressure without actually solving numerically the Poisson equation. Based on an inverse kinetic 
approach recently proposed for the incompressible Navier-Stokes equations we intend to prove that 
an exact evolution equation can be obtained which advances in time self-consistently the fluid pres- 
sure. The new equation is susceptible of numerical implementation in Lagrangian CFD simulation 
codes. 
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PACS numbers: 47.27.Ak, 47.27.Eq, 47.27.Jv 

An important aspect of computational fluid dynam- 
ics is related to the determination of the fluid pressure 
in isothermal incompressible fluids and particularly to 
the construction of an exact evolution equation for the 
fluid pressure which replaces the Poisson equation. This 
amounts to transform an elliptic type fluid equation into 
a suitable hyperbolic equation, a result which usually is 
reached only by means of an asymptotic formulation. In 
this paper we intend to show that an exact solution to 
this problem is possible when the evolution of the fluid 
fields is described by means of a suitable dynamical sys- 
tem, to be identified with the so-called Navier-Stokes 
(N-S) dynamical system Besides being a still un- 
solved mathematical problem, the issue is relevant for at 
least two reasons: a) the proliferation of numerical algo- 
rithms in computational fluid dynamics which reproduce 
the behavior of incompressible fluids only in an asymp- 
totic sense (see below) ; b) the possible verification of con- 
jectures involving the validity of appropriate equations 
of state for the fluid pressure. Another possible motiva- 
tion is, of course, the ongoing quest for efficient numer- 
ical solution methods to be applied for the construction 
of the fluid fields {/?, V,p} , solutions of the initial and 
boundary- value problem associated to the incompressible 
N-S equations (INSE). For dcfinitcncss, it is convenient 
to recall that this is defined by the continuity, N-S and 
isochoricity equations 

JVV = 0, (1) 
V • V = 0, (2) 

where the mass density p and the fluid pressure p are 
required to satisfy the inequalities 



p(r,t) > 0, 
P( T J) = Po > 0, 



(3) 
(4) 



a suitably smooth volume force density acting on the fluid 
element, p a is the constant mass density and fi > is the 
constant fluid viscosity. Equations (JT|)- (J4j are assumed 
to admit strong solutions in an open set f2 x J, with 
f! C I 3 the configurations space (defined as the subset 
of M 3 where p(r,t) > 0) and I C K a possibly bounded 
time interval. By assumption {p, V,p} are continuous in 
the closure f2. Hence if in fl X I, f is at least C( 1,0 '(f2 x 
/), it follows necessarily that {p, V,p} must be at least 
C (2 ^(n x I), while the fluid pressure and velocity must 
satisfy respectively the Poisson and energy equations 



V 2 p = -V -f-pV-(V- VV), 



V-NV = o. 



(•5) 



(6) 



Here N is the N-S operator iVV =p^V+Vp+{-fiV 2 V, 



with -jjj = Jj+V • V the convective derivative, f denotes 



It is well known that the choice of the Poisson solver 
results important in numerical simulations, since its effi- 
cient numerical solution depends critically on the num- 
ber of modes or mesh points used for its discretization 
(see also Ref.Q and references therein indicated). In 
turbulent flows this number can become so large to ef- 
fectively limit the size of direct numerical simulations 
(DNS) 0| . This phenomenon may be worsened by the 
algorithmic complexity of the numerical solution meth- 
ods adopted for the Poisson equation. For this reason 
previously several alternative approaches have been de- 
vised which permit to advance in time the fluid pressure 
without actually solving numerically the Poisson equa- 
tion. Some of these methods are asymptotic, i.e., to ad- 
vance in time the fluid pressure they replace the exact 
Poisson equation with suitable algorithms or equations 
which hold only in an asymptotic sense (neglecting suit- 
ably small corrections), others are exact solvers, i.e., pro- 
vide in principle rigorous solutions of INSE (and Poisson 
equation) . The first category includes the pressure-based 
method (PBM) 0, the Chorin artificial compressibility 
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method (ACM) [6j, the so-called preconditioning tech- 
niques [? ], all based on ACM, and kinetic approaches, 
of which a notable example is provided by the so-called 
Lattice-Boltzmann (L-B) methods (for a review see for 
example Ref. 11[ and references therein indicated) . PBM 
is an iterative approach and one of the most widely used 
for incompressible flows. Its basic idea is to formulate 
a Poisson equation for pressure corrections, and then 
to update the pressure and velocity fields until the iso- 
choricity condition ([2|) is satisfied in a suitable asymp- 
totic sense. The ACM approach and the related precon- 
ditioning techniques, instead, are obtained by replacing 
the Poisson and N-S equations with suitable parameter- 
dependent evolution equations, assuming that the fluid 
fields depend on a fictitious pseudo-time variable r. In 
dimensionless form the evolution equation for the pres- 
sure becomes in such a case £ 2 -§^P + V • V =0, where 
e 2 > is an infinitesimal parameter. Manifestly this 
equation recovers only asymptotically, i.e., for e 2 — > 0, 
the exact isochoricity condition j2|) . Introducing the fast 
variable r = r/e 2 , this implies that the fluid fields must 
be of the form V(r, t,r),p(r, t,r) and should be assumed 
suitable smooth functions of r. Therefore, for prescribed 
finite values of e 2 ( to be assumed suitably small), this 
equation permits to obtain also an asymptotic estimate 
for the fluid pressure p(r,t). This is expressed by the 
equation 

p(r,t) = Jim p(r,t,r) ^ (7) 

r — >oo 

<* p(r,t,T = 0)-f dr'V-V(r,i,r'), 
o 

where r* >> 1 is suitably defined and p(r,t,T = 0) de- 
notes some initial estimate for the fluid pressure. Several 
implementations on the Chorin algorithm are known in 
the literature (see for example Refs.@, IE EJ E3]). Cus- 
tomary L-B methods are asymptotic too since they re- 
cover INSE only in an approximate sense; moreover typ- 
ically they rely on the introduction of an equation of 
state for the fluid pressure, for example, the equation 
of state of an ideal gas, or more refined models based 
on so-called non- ideal fluids This assumption, how- 
ever, generally requires that the effective Mach-number 
characterizing the L-B approach, defined by the ratio 
M e ff = V sup /c (with c denoting the discretized veloc- 
ity of the test particles and V sup the sup of the velocity 
field at time f), must result suitably small. As a con- 
sequence, in typical L-B approaches the fluid pressure 
can only be estimated asymptotically. However, there 
are other numerical approaches which in principle pro- 
vide exact Poisson solvers. These include the so-called 
spectral methods in which the fluid fields are expanded 
in terms of suitable basis functions. Significant exam- 
ples are the pure spectral Galerkin and Fourier methods 
[1 3l | as well as the nonlinear Galerkin method [14| , which 



are typically adopted for large-scale turbulence simula- 
tions. In these methods the construction of solution of 
the Poisson equation is obtained analytically. However, 
the series-representation of the fluid fields makes diffi- 
cult the investigation of the qualitative properties of the 
solutions, such - for example - the search of a possible 
equation of state or an evolution equation for the fluid 
pressure. 

Another approach which provides in principle an ex- 
act Poisson solver is the one recently proposed by Ellcro 
and Tessarotto 0, Q| , based on an inverse kinetic theory 
for INSE. This approach, recently applied also to quan- 
tum hydrodynamic equations [151 ]. permits to represent 
the fluid fields as moments of a suitably smooth kinetic 
distribution function /(x, i) which obeys an appropriate 
inverse Vlasov-type kinetic equation: 



d d 



(8) 



Here X(x,i) = {v,F} and x = (r, v) er C Q x M 3 is 
the state vector generated by the vector field X, v is the 
kinetic velocity, while F(x,£) is an appropriate mean-field 



17], it has been 



force obtained in Ref. In Refs. [If 
proven that F(x,i) can be uniquely prescribed, in partic- 
ular, in such a way that: 

• All the fluid equations are obtained from appropri- 
ate moments of Eq.©. As a consequence, the fluid 
equations as well as the initial and boundary con- 
ditions for the fluid fields are satisfied identically. 

• The time evolution of the kinetic distribution func- 
tion, T t)to /(x ) = /(x(i),t), is determined by the 
classical dynamical system associated to the vector 
field X, i.e., 



5 x = X(x,t) 
x(t D ) = x 



(9) 



(N-S dynamical system) which must hold for arbi- 
trary initial conditions x Q = (r ,v ) G T. 

The solution of @, x(i) = T tjto x , which de- 
fines the N-S evolution operator T t . to , determines 
uniquely a set of curves {x(i)} = {x(i), Vt G J} x 
obtained for arbitrary (x ,t ) G T x /, which can 
be interpreted as phase-space Lagrangian trajec- 
tories associated to a set of fictitious "test" par- 
ticles. Their projections onto the configuration 
space, denoted as configuration- space Lagrangian 
trajectories, are defined by the curves {r(i)} = 
{r(t) = T tjto r ,Vf G /} . By varying their initial 
conditions, in particular r G fi, the curves {r(i)} 
can span, by continuity, the whole set f2. 

The fluid pressure p(r, t) is defined by 

p(r,t)= Pl (r,t)-p (t), (10) 
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(to be regarded as a constitutive equation for 
p(r, i)), where pi(r, t) is the kinetic pressure 
Pi(r,t) — J dv^u 2 f(x,t), p is denoted as re- 
duced pressure, while u is the relative velocity 
u=v - V(r,*). 

• By definition, the reduced pressure p is solely a 
function of time, to be assumed suitably smooth 
and prescribed. Both p a (t) andpi(r, t) are strictly 
positive, while p (t) in fl x I is subject to the con- 
straint pi(r, t) — p {t) > 0. 

• A particular solution of the inverse kinetic equation 
([5]) is provided by the local Maxwellian distribution 
/M(x,t;V, Pl ) = _^_exp{-r 2 } [where Y 2 = 



and vf, 

vth z tn 

field F reads: 
F(r,v,i) 



2pi/p ]. In such a case, the vector 



1 



Vp{ 
P \.Pi 



a--N V+^A oPl - 
P 2 

3 
~ 2 



(11) 



where a denotes the convective term a=^u ■ 
VV+iW • u, £ is the relative kinetic energy den- 
sity £=pu 2 /2, while Nq and Aq are the differ- 
ential operators -/V V = — f(r, V,t) + pV 2 V and 
A oPl (r,t) = i Pl + V • (Vpi)J. For an arbi- 
trary and suitably smooth distribution function 
/(x,i), the form of the vector field F satisfying these 
hypotheses has been given in Refs. [3, Eg. 

An interesting issue is related to the consequences of 
the constitutive equation (fTOf and of the N-S dynamical 
system generated by the initial value-problem ([9]). In 
this Letter we intend to prove that the fluid pressure 
p(r, t) obeys an exact partial-differential equation which 
uniquely determines is time evolution. This is obtained 
by evaluating its Lagrangian derivative along an arbi- 
trary configuration-space Lagrangian trajectory {v(t)} 
generated by the N-S dynamical system. The result can 
be stated as follows. 

Assuming that the initial-boundary value problem asso- 
ciated to INSE admits a suitably strong solution {p, V,p} 
in the set x I, the following statements hold: 

A) If x(t) is a particular solution of Eq. (0|) which 
holds for arbitrary r(t) 6 f2 and t € I, along each phase- 
space Lagrangian trajectory {x(t)} defined by Eq. ([9j) 
the scalar field £(r,<) = £/ P i obeys the exact evolution 
equation 



d , 1 

which holds for arbitrary initial conditions x D 



(12) 



and £ Q 



2 Pl (r ,t ) ' 



with u = v — V(r , t Q ). Here is 



d_ 

di 



while all quantities (u, E and p\) are evaluated along an 
arbitrary phase- space trajectory {x(i)}. 

B) Vice versa, if the solutions x(f)=(r(f), v(f)) and 
£(i) of Eqs.^j, are known for arbitrary initial con- 
ditions x = (r o ,v D ), u D = v G - V(r ,i ) and £ D = 

2 

2p! (r " 1 ) ) an d for all (r,i) E fix I, it follows necessar- 
ily that in fl x I, {p, V,p} satisfy identically INSE. 
PROOF 

Let us first prove statement A) , namely that INSE and 
the N-S dynamical system imply necessarily the validity 
of Eq. (fT2|) . For this purpose we first notice that by con- 
struction Eq.© admits a unique solution x(t) for arbi- 
trary initial conditions x„ = (r Q , v Q ) £ T. while the same 
equation can also be expressed in terms of the relative 
velocity u= v — V(r,i). This yields 



d u = F- W(r ' t) -u.VV(r,t) (13) 



dt 



Dt 



Upon invoking the N-S equation Jl} and by taking the 
scalar product of Eq. (|13p by pu, this equation implies 



dt 



£=u-Vp 



{- 


4} 


£ 


I Pi 




H 

Pi 



Pi + V • (Vp) 



(14) 



which gives 



2pi 



u-Vp+fV-V. (15) 



the Lagrangian derivative 4^ = J^+v-V + F-^, £(r, t), 



As a consequence of the isochoricity condition (J2J this 
equation reduces identically (i.e., for arbitrary initial con- 
ditions for the dynamical system) to Eq. (fT2l . B) Vice 
versa, let us assume that the solutions x(£) = (r(t), v(t)) 
and of Eqs.©, (fl"2")) are known for arbitrary initial 

2 

conditions x Q G T and £ G = 2 p 1 ^r ° t ) • l n this case ^ 
follows the fluid fields necessarily must satisfy INSE in 
the whole set ft X I. It suffices, in fact, to notice that 
by assumption the evolution operator T tjto is known. 
This permits to determine uniquely the kinetic distri- 
bution function at time t, which reads [l[ /(x(i),i) = 
f(x ,t )/J(x(t),t), where J(x(t),t) is the Jacobian of 
the flow x Q — > x(t). Hence, also its moments are uniquely 
prescribed, including both V(r,t) and p(r,t), in such a 
way that they result at least C' 2 ' 1 ^^ x 7). The inverse 
kinetic equation (jSJ) , thanks to the special form of F as 
given by Eq. ifTTj) ensures that the N-S equation is sat- 
isfied identically in SI x I [l|. Moreover, since Eqs. (fT3|) 
and (|15p are by assumption fulfilled simultaneously, it 
follows that both the isochoricity condition ([2]) and the 
Poisson equation [Eq.©] must be satisfied too in £1 x I. 
This completes the proof. 

As a basic implication, if the fluid velocity is assumed 
to satisfy both the N-S equation and isochoricity con- 
dition, the mass density satisfies the incompressibility 
condition Q, while {x(i)} is an arbitrary trajectory of 
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the N-S dynamical system, it follows that Eq. (fT2|) deter- 
mines uniquely the time-advancement of the fluid pres- 
sure. Hence, it provides an evolution equation for the 
fluid pressure, which by definition is equivalent simul- 
taneously to the isochoricity condition and to the Pois- 
son equation. This equation can in principle be used 
to determine the fluid pressure at an arbitrary position 
r £ f2. However, since any given position can be reached 
by infinite phase-space (and also configuration-space) La- 
grangian trajectories, it is sufficient to sample the config- 
uration space by a suitable subset of Lagrangian trajec- 
tories (test particles), obtained by prescribing the initial 
condition x„. 

The physical interpretation of the pressure evolution 
equation is elementary: it yields an unique prescription 
for the Lagrangian time derivative of the fluid pressure, 
which is defined in the frame which is locally co-moving 
with a test particle of state x(t) and velocity v(£).In par- 
ticular, it is obvious that the specification of the initial 
kinetic velocity v Q remains essentially arbitrary, as well 
as the definition of the reduced pressure p (t). This 
means that both the dimcnsionless ratios My = V/ |v D | 
and M p = p/p a , to be denoted as velocity and pressure ef- 
fective Mach numbers, remain essentially arbitrary. As 
a consequence it is possible, in principle, to construct 
asymptotic solutions of Eq. ((T2|) based on low effective- 
Mach numbers expansions, i.e., for which My, M p <C 1. 
As an illustration, let us prove that an approximate so- 
lution of this type can be obtained for p\ (and hence 
for p) in the so-called diffusive approximation, i.e., by 
considering a subset of velocity space in which by as- 
sumption at time t the initial relative velocity |u| t 
and the relative kinetic energy £ =p a u 2 /2 are assumed 
to satisfy the Mach-number orderings My ~ 0(5) and 

1 

^ 5> 



O(d ), being 6 < 1. 



These imply £ 



l u lt ~ jt/j. It follows that in an infinitesimal time 
interval [t a , t\ = t + At], assuming At ~ 0(8), there 



results gjU 



P ^Pl 



(diffusive approxima- 



l + 0(S 1/2 ) 

tion) which yields, by integrating it in the Euler approx- 
imation, 



u(t) - u(t ) S -Vp—At \l + 0(5 1/2 ) 
P Pi L 



(16) 



In the same approximation the relative kinetic energy at 
time t becomes 



£(t) 



£(to) 



1 - u(t„) -V In pi (t„) At 



1 + 0{S 1/2 ) . (17) 



As a consequence, Eq. (fT2)) can now be used to advance 
in time p\. In fact, integrating it and invoking again the 
Euler approximation, yields 



£_ 

Pi 



£_ 

Pi 



s --u(f ).Vlnpi(f )At, (18) 



which delivers an equation for pi (t) . We stress that these 
features are potentially important for the construction of 
possible numerical algorithms based on Eq. (fT2|) . There- 
fore, the pressure evolution equation can in principle be 
adopted for the development of Lagrangian particle sim- 
ulation methods in fluid dynamics. These developments 
will be the object of future investigations. 
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